from brian import *
from brian.library.random_processes import *
from brian.library.synapses import  *
import SuperNeuron as SuperNeuronCl
from MyConstants import *



class InNeuronFactory (SuperNeuronCl.SuperNeuron):
    eqs=Equations('''
    dh/dt= (hinf-h)/Tauh : 1
    dn/dt= (ninf-n)/Taun : 1
    dmnap/dt = (mnapinf-mnap)/Taumnap : 1
    dhnap/dt = (hnapinf-hnap)/Tauhnap : 1
    dm/dt = (minf-m)/Taum : 1
    dvm/dt = 1/C*(-INa-INap-IK-Ileak-ISynI-ISynE-Irand) : mV
    
    Tauh = 30/(exp((vm/mV+50)/15)+exp(-(vm/mV+50)/16))*ms : ms
    Taun = 7/(exp((vm/mV+40)/40)+exp(-(vm/mV+40)/50))*ms : ms
    Tauhnap = 1200/(cosh((vm/mV + 59)/16))*ms :ms
    
    minf = 1/(1+exp((vm-Vhm)/Sm)) : 1
    hinf = 1/(1+exp((vm-Vhh)/Sh)) : 1
    ninf = 1/(1+exp((vm-Vhn)/Sn)) : 1
    mnapinf = 1/(1+exp((vm-Vhmnap)/Smnap)) : 1
    hnapinf = 1/(1+exp((vm-Vhhnap)/Shnap)) : 1
    
    INap = GNap*mnap*hnap*(vm-ENa) : mA*umetre**-2
    INa = GNa*m**3*h*(vm-ENa) : mA*umetre**-2
    IK = GK*n**4*(vm-EK) : mA*umetre**-2
    Ileak = gleak*(vm-Eleak) : mA*umetre**-2
    
    gMlrE = gEd * drive : siemens*metre**-2
    
    ISynI = (gMlrE+ipsp)*(vm-Ei) : amp*metre**-2
    ISynE = (gMlrE+epsp)*(vm-Ee) : amp*metre**-2
    
    GNa :  msiemens*cm**-2
    GNap :  msiemens*cm**-2
    GK :  msiemens*cm**-2
    gleak :  msiemens*cm**-2
    Eleak : mV
    
    drive : 1
    ''')
    def init(self, N_RG=0,N_PF=0,N_IN=0,N_INRG=0,synapse='exp',Evariable=1,rand=False):
        self.add_synapses(synapse,tau1=tausyn)
        self.add_rand_curr(rand)
        self.rand=rand
        self.Evariable = Evariable
        self.N_RG = N_RG
        self.N_PF = N_PF
        self.N_IN = N_IN
        self.N_INRG = N_INRG
        self.Group=NeuronGroup(self.N_RG+self.N_PF+self.N_IN+self.N_INRG,model=self.eqs,
            threshold=EmpiricalThreshold(threshold=0*mV,refractory=5*ms),
            implicit=True)
        self.SG_RG = self.Group.subgroup(self.N_RG)
        self.SG_PF = self.Group.subgroup(self.N_PF)
        self.SG_IN = self.Group.subgroup(self.N_IN)
        self.SG_INRG = self.Group.subgroup(self.N_INRG)
        #print self.eqs
        self.reinit()
        return self.Group
    
    def reinit(self):
        if self.N_RG != 0:
            
            self.SG_RG.GNa = 30 * msiemens*cm**-2
            self.SG_RG.GNap = 0.25 * msiemens*cm**-2
            self.SG_RG.GK = 1 * msiemens*cm**-2
            self.SG_RG.gleak = 0.1 * msiemens*cm**-2
            self.SG_RG.Eleak = (-64+randn(self.N_RG)*0.64*self.Evariable) * mV
            self.SG_RG.drive = 0.5
            
        
        if self.N_PF != 0:
            
            self.SG_PF.GNa = 30 * msiemens*cm**-2
            self.SG_PF.GNap = 0.1 * msiemens*cm**-2
            self.SG_PF.GK = 1.2 * msiemens*cm**-2
            self.SG_PF.gleak = 0.1 * msiemens*cm**-2
            self.SG_PF.Eleak = (-64+randn(self.N_PF)*0.64*self.Evariable) * mV
            self.SG_PF.drive = 0.5
            
        if self.N_IN != 0:
            
            self.SG_IN.GNa = 120 * msiemens*cm**-2
            self.SG_IN.GNap = 0 * msiemens*cm**-2
            self.SG_IN.GK = 100 * msiemens*cm**-2
            self.SG_IN.gleak = 0.51 * msiemens*cm**-2
            self.SG_IN.Eleak = (-64+randn(self.N_IN)*3.2*self.Evariable) * mV
            self.SG_IN.drive = 0
            
        if self.N_INRG != 0:
            
            self.SG_INRG.GNa = 120 * msiemens*cm**-2
            self.SG_INRG.GNap = 0 * msiemens*cm**-2
            self.SG_INRG.GK = 100 * msiemens*cm**-2
            self.SG_INRG.gleak = 0.51 * msiemens*cm**-2
            self.SG_INRG.Eleak = (-57.5+randn(self.N_INRG)*2.875*self.Evariable) * mV
            self.SG_INRG.drive = 0
            
         
        self.Group.vm = -70*mV
        self.Group.mnap = 1/(1+exp((self.Group.vm-Vhmnap)/Smnap))
        self.Group.hnap = 1/(1+exp((self.Group.vm-Vhhnap)/Shnap))
        self.Group.m = 1/(1+exp((self.Group.vm-Vhm)/Sm))
        self.Group.h = 1/(1+exp((self.Group.vm-Vhh)/Sh))
        self.Group.n = 1/(1+exp((self.Group.vm-Vhn)/Sn))
            
